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Multiple time scale molecular dynamics enhances computational efficiency by updating slow mo- 
tions less frequently than fast motions. However, in practice the largest outer time step possible 
is limited not by the physical forces but by resonances between the fast and slow modes. In this 
paper we show that this problem can be alleviated by using a simple colored noise thermostatting 
scheme which selectively targets the high frequency modes in the system. For two sample problems, 
flexible water and solvated alanine dipeptide, we demonstrate that this allows the use of large outer 
time steps while still obtaining accurate sampling and minimizing the perturbation of the dynamics. 
Furthermore, this approach is shown to be comparable to constraining fast motions, thus providing 
an alternative to molecular dynamics with constraints. 



I. INTRODUCTION 

Over the past two decades atomistic simulation of 
chemical, material and biological systems has become a 
routine and important tool for understanding, analyzing 
and predicting experiment. Time scales of interest now 
often span into the microsecond range so as to monitor 
processes such as protein folding and transport through 
membranes. However, the shortest atomic time scales 
in these systems, typically for bonded interactions such 
as stretches and bends, are of the scale of femtoseconds 
and therefore time steps of this order are required to sta- 
bly evolve the system. As a result, billions of steps are 
needed to reach the time scales of interest, making such 
simulations a formidable challenge. 

Chemical systems typically involve interactions occur- 
ring on many time scales ranging from rapidly varying, 
but cheap to calculate, bonded interactions to slow, but 
expensive, long range electrostatics. Multiple time scale 
methods can be used to exploit this separation in time 
scales by updating the slow interactions less frequently 
than the fast interactions, in principle allowing signifi- 
cant computational savings to be achievedPSl However, 
the maximum outer time step which can be obtained is 
in practice lim ited by the resonance between the slow 
and fast modesP^ This results in energy building up in 
the high frequency modes, raising the temperature of the 
system and leading to unstable trajectories and incorrect 
sampling. As a result the fastest mode in the system 
still dictates how frequently the slowest interactions must 
be calcul ated , thereby limiting the maximum obtainable 
speed-upP^ 

A commonly used approach to delay the resonance bar- 
rier is to remove the high frequencies in the system by 
constraining themf^El For example in water one can con- 
strain the oxygen-hydrogen bond length and intramolec- 
ular hydrogen-hydrogen distance. This allows a larger- 
outer time step to be used but requires prior knowledge of 
the coordinates that comprise the fast modes. In simple 



systems these coordinates may be known, but in complex 
biological and materials systems this is not always the 
case and one could risk constraining a degree of freedom 
which has vital importance to mechanism or function. 
The MOLLY method can be viewed in a similar fashion 
since it requires prior specification of coordinates of fast 
modes which are then used to filter o ut the destabilizing 
components of the slow forces EHfiUil 

It is now well established that large time steps can 
be achieved while retaining full flexibility by coupling 
the system to a bathPSHHI We note that the bath serves 
two purposes: to remove energy which builds up in the 
modes and to disrupt the high frequency modes from res- 
onating with the lower frequencies. The most common 
choice is to couple each atom in the system to a white 
noise Langevin ( WNL) bath which acts uniformly across 
the spectrum of the systemP^Hll However, as has been 
pointed out previously,^ the strength of the system-bath 
coupling needed to stabilize large time steps significantly 
disrupts the motion of the slow modes thereby greatly 
hindering diffusive and orientational motion. As we will 
show this leads to a situation where any computational 
speed-ups gained by increasing the outer time step are 
largely outweighed by the decrease in the rate at which 
different configurations are explored. To avoid this one 
would ideally like to couple strongly to high frequency 
motion in the system while leaving low frequencies un- 
perturbed. This can be achieved by transforming to the 
normal modes of the system at each time step and then 
evolving using a strong coupling to high frequencies with 
weak coupling to low ones. Indeed, it has been shown for 
a simple biological system in implicit solvent that this 
approach is very successful at removing resonance issues 
while preserving dynamics.^ However, in practice per- 
forming normal mode transformations at each step for 
large systems is prohibitively expensive. 

In this paper we attempt to reconcile the simplicity of 
the WNL approach with the effectiveness of the targeted 
normal mode approach by using tailored colored noise. 
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Recent work has made this possible by showing how 
the generalized Langevin equation (GLE) can be used to 
thermostat molecular dynamics simulations using an ex- 
tended WNL formalismEsmni Unlike white noise, colored 
noise can be tailored to have a frequency dependent cou- 
pling which allows for much greater flexibility in its ap- 
plication. For example a simple colored noise thermostat 
was designed for use in Car-Parrinello ab-initio molecular 
dynamics simulations which only targets atomic motion 
while allowing the high frequency fictitious electronic de- 
grees of freedom to evolve freelyP^ In this work we will 
demonstrate how colored noise can instead be used to de- 
sign a bath which heavily damps high frequencies while 
leaving low frequencies largely unaffected. The simple 
scheme which results allows the resonance barrier to be 
postponed facilitating the use of large outer time steps 
while yielding accurate sampling and minimal impact on 
the dynamics. Applications of this approach to flexi- 
ble water and an aqueous solution of alanine dipeptide 
demonstrate that significant increases in computational 
efficiency can be achieved while requiring little a priori 
knowledge of the system. 

The outline of the paper is as follows. Section II briefly 
reviews the colored noise thermostatting approach and 
shows how colored noise profiles can be constructed in 
a transparent way by combining simple analytic forms. 
The implementation of the GLE thermostat in a stan- 
dard reference system propagator algorithm (RESPA) 4 
is then discussed. Section III outlines the simulations 
performed using this scheme to calculate the static and 
dynamic properties of a fully flexible model of liquid wa- 
ter and an aqueous solution of alanine dipeptide. Section 
IV discusses the results of these applications and Sec. V 
concludes. 



In the case where the random force is uncorrelated, 
(R(t)R(0)) — b 2 5(t), the white noise Langevin equation, 



x = p/m 

p = f(x)-jp + bt;(t) 



(4) 
(5) 



is recovered. Here £(t) is a Gaussian Markov process with 
unit variance and the fluctuation dissipation theorem in 
Eq. ([3| reduces to, 



2mk B T^ = b 2 



(6) 



which is commonly used as a tool to thermostat molecu- 
lar dynamics simulations. 

A colored noise thermostat can be implemented by ex- 
actly mapping Equation [2] onto a Markovian dynamics 
in an extended space consisting of a set of auxiliary mo- 
menta, s, that are coupled to the syst em mom entum, p, 
in the presence of a Markovian bathPtHTMH] 

The equa- 
tions of motion are given by: 



p/m 

(V 



■B£(t) 



(7) 
(8) 



where £(t) is a vector of uncorrelated Gaussian noise. 
The drift (or damping) matrix (T), may be related to the 
diffusion matrix (B) by a recasted fluctuation-dissipation 
theoremPEl 



mk B T T 



BB 



The matrix elements of T, 



p _ / Ipp Ips 

^fsp T ss 



(9) 



(10) 



II. THEORY 



Colored noise thermostats 



For a particle of mass m with position x and momen- 
tum p moving on a one dimensional potential energy sur- 
face V(x), the generalized Langevin equation ispHES 



p/m 



J drK(t - t)p(t) + R(t) 



(1) 
(2) 



determine the form of the memory kernel that acts on 
the system via 



K(t) = 2 lpp 5(t) 



ips 



-I*|r 3 



Is 



(11) 



If the extended momenta are uncoupled from the system 
equation, Eq. (8) reduces to the white noise Langevin 
equation (Eq. (5)) with j vv — 7. In the next section 
we will describehow Eq. ( |lT| ) can be used to design a 
drift matrix that yields a colored noise profile such that 
the high frequency modes are strongly coupled and slow 
modes weakly coupled to the bath. 



where K(t) is the memory kernel, R(t) is a non- 
Markovian "colored" random force and f(x) = 
— dV(x)/dx is the force due the potential. The 
fluctuation-dissipation theorem dictates that for an equi- 
librium system at temperature T, K (t) and R(t) are re- 
lated by, 



mk B T K{t) = (R(t)R(0)) 
where k B is the Boltzmann constant. 



(3) 



B. Designing the colored noise profile 

The effect of the bath on the system is determined 
by the details of K[t). Using the extended white noise 
Langevin framework outlined above the effect of the bath 
and the form of the memory kernel is determined by cus- 
tomizing the drift matrix T (see Eq. (11)). 



The expression given in the right hand side of Eq. ( 11 ) 



consists of two terms, a simple white noise component of 
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strength j pp and one that depends upon the coupling of 
the auxiliary momenta. This second term is a scalar and 
is therefore invariant to the choice of basis. The time de- 



pendence of the second term of Eq. (11) is contained in 
the matrix exponential of submatrix, T ss and it is there- 
fore transparently expressed in its eigenbasis. The eigen- 
decomposition of the submatrix yields T ss — UAU -1 , 
where = XjSij is the diagonal eigenvalue matrix. 
Matrix U is comprised of the set of eigenvectors in its 
columns and may be used to rotate, T ss into the eigen- 
basis. The memory function is now expressed as, 



K{t) = 2 lpp 5(t) 



ItlAjj- 



Isp- 



(12) 



In order to recover correct equilibrium behavior in the 
t — > oo limit, it is a necessary but not sufficient condition 
that the re al par t of the eigenvalues of T ss are chosen to 
be positive? 18 * 23 ! The resultant kernel consists of a linear 
combination of exponential decaying forms, arising from 
the real eigenvalues, and damped oscillatory ones from 
the complex conjugate pairs. 

In considering the effect of a colored noise bath it is 
physically intuitive to consider the memory spectrum, 
K(uj), which is the Fourier transform of the memory ker- 
nel. Given that the memory function may be expressed 
as a linear combination of a white noise term and a set of 
exponential forms with real or complex conjugate pairs of 
eigenvalues, the memory spectrum is readily computable. 
For a given eigenvalue Xj = a,j +ibj , the Fourier transform 
of the corresponding exponential function is a Lorentzian, 



dt 



-iu)t^ibjt^- 



b 3 y 



(13) 



The width and center of the Lorentzian is given by the 
real and imaginary parts of the corresponding eigenvalue, 
respectively. Therefore, the memory spectrum is express- 
ible as a constant 2^ pp from which a set of Lorentzian 



functions are added or subtracted. In order to gener- 
ate the correct equilibrium distribution, forms m ust be 
chosen such that K{ui) is always greater than zeroP^l 
Although the form of the memory function as given 
by Eq. \12\ appears simple, "translating" a chosen form 
into a stable drift matrix is non-trivial. This is due to the 
fact that although the basis in which the submatrix T ss 
is expressed is arbitrary with respect to the form of the 
memory function, this choice of basis is crucial for the 
generation of stable dynamics within the colored noise 
thermostat implementation. Namely, a requirement of 
the integration of the equations of motion is that the 
symmetric part of T must be positive definite (see Sec. 



II C). Indeed, it is possible to choose forms of the kernel 
based on Eq. ( 12 ) such that K(u>) > while violating 



this condition. 

Despite these difficulties, it is possible to construct 
drift matrices which correspond to a memory spectrum 
of a desired form where the matrix elements are trans- 
parently relatable to the the time dependence, namely 
the eigenvalues of T ss are explicit parameters of T. Drift 
mat rices of this sort have already appeared in the litera- 
£ ure LL6i!9] an( j their details are summarized in Appendix 
[A) We now present such a form that is tailored for the 
problem at hand. 

In selecting the details of the colored noise profile it is 
instructive to first consider the memory spectrum, which 
provides information about how strongly the modes of 
the system exchange energy with the bath at a given fre- 
quency.^ Additionally, the oj — value of the memory 
spectrum is inversely proportional to the diffusion con- 
stant of a free particle (/(x) = 0) coupled to the bath. 
Therefore, if a profile is required that couples the bath 
strongly to high frequency modes and more weakly to 
low ones, one can start by constructing a memory kernel 
whose spectrum is small at low frequencies and large near 
the frequencies where strong coupling is desired. The fol- 
lowing drift matrix, 



7oo 3 1/4 v /w(7 00 - 70 ) ij a/£(7oo - To) 



3^ ^(700 - 7o~) "VS Co ) (14) 

-3173-^(700 - To) 

corresponds to a memory spectrum which possesses such a form. The associated memory kernel and spectrum are 
given by the following equations: 

K(t) = 2 7 ^(()-2-^( 7 „-7„)e--^l'l/ 2 cosM/2) (15) 

A ' M " ^ - 2 ^3 h '~ - ( to-n-f-f)' + W^W ) ■ <16) 



The memory spectrum of Eq. ( 16 ) is shown schemati- cally in Figure [T] It can be understood as a white noise 
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of strength 7^ from which a pair of Lorentzians of width 
o)v / 3/2 centered at uj — ±uj/2 is subtracted. In this way 
the Lortenzian terms can be thought of as partially "can- 
celing out" the white noise component at low frequencies. 
The shape is such that its value is 270 at uj = and ap- 
proaches a value 2700 as uj — > 00. In order to fulfill our 
design requirements of the noise profile, we take the pa- 
rameter 70 to be an arbitrarily small (non-zero) number. 
The parameter u> is related to the frequency at which the 
memory spectrum begins to rise rapidly. When the pa- 
rameters are chosen to have positive values, K{uj) is an 
even, positive definite function, as is required to generate 
the correct equilibrium distribution. 

The colored noise profile expressed in Eqs. (14), (151, 



and ( 16 1 provides a transparent and readily tunable form 
for the drift matrix of the colored noise thermostat with 
an associated memory spectrum of the shape that we 
desire. However the coupling of the bath to system is 
non-local, and although intuition can inform our choice 
of the shape of the memory spectrum, one cannot a pri- 
ori determine the exact values of the parameters 700 
and uj which adequately couple to high frequency modes 
and simultaneously minimize the perturbation on low fre- 
quency modes. 

In previous work, colored noise thermostats have been 
successfully tuned according to the energy relaxation 
times of harmonic oscillators coupled to the thermo- 
stat.^ However in this work we find it more appropri- 
ate to estimate the disturbance of the spectrum of a 
test harmonic oscillator engendered by the colored noise 
thermostat. The velocity autocorrelation function (and 
hence spectrum) for an oscillator with frequency ujq in 
the presence of the colored noise t herm ostat may be ex- 
actly computed by matrix algebra? 18 ! 27 ! The spectrum is 
both broadened and shifted with respect to the free os- 
cillator result, and the size of these differences provides 
some measure of the strength of the system-bath cou- 
pling. For frequencies that are small compared to the 
bath parameter uj and when 70 is negligible, the ratio of 
the peak position of the system-bath spectrum (wp) to 
that of the free spectrum (ujq) may be estimated by the 
following expression (see Appendix |B| , 



UJp 

w 



1 



O 



1 



dV3 



(17) 



It is important to note that this estimate depends on 
the value of the friction in the high frequency limit, and 
is a function of the ratio 7^ jUj. This underscores the 
non-locality in frequency space of the system-bath inter- 
action. 

Equation |17| provides an estimate of the impact of the 
colored noise thermostat on the low frequency modes. 
In order to ensure that the high frequency modes are 
sufficiently damped, the parameter 700 is set to be large. 
A lower bound is provided by the white noise friction 
that is required to stabilize any resonance instabilities 
that are present in the system. Fine tuning can then 



matrix no. 


7oo (ps L ) 


70 (ps L ) 


& (ps L ) 


Up/ UJQ 


GLE - 12fs 


83.33 


0.01 


300.0 


0.93 


GLE - 16fs 


125.0 


0.01 


100.0 


0.76 


GLE - 20fs 


200.0 


0.01 


75.0 


0.63 



TABLE I: The parameters sets of the drift matrix (Eq. ( 14 ) ) 
which are utilized in this study. Parameter sets are labeled 
according to the outer time step which they are used in con- 
junction with (See Sec. 111). 



be accomplished by testing the performance of a set of 
noise profiles on a realistic system in conjunction with 
a multiple time scale integrator. However, since most 
biomolecular systems have similar spectral features the 
profiles presented here should provide good performance 
in a broad range of studies without reparameterization. 
The parameters for the colored noise profiles that will be 
utilized in this study are given in Table [I] 

In order to gauge the impact of the colored noise ther- 
mostat, the spectrum of a flexible water modePSl i n the 
presence of a bath as defined by parameter set GLE - 12 fs 
in Table [T] is plotted in Figure [2j The result is compared 
to microcanonical and white noise Langevin dynamics 
with a friction that is equivalent to the uj — > 00 limit 
of the colored profile. It can be readily seen that, when 
compared with the microcanonical dynamics, the lower 
frequency modes are less disturbed than those related to 
bending and stretching. It is the intramolecular modes 
that have been targeted for damping and one can see that 
the impact of the colored bath on the oxygen-hydrogen 
stretch at sa 680 ps -1 is comparable to that engendered 
by the high-friction white noise bath. Furthermore we 
note that the ratio of the low frequency peak positions 
of the thermostatted spectra to the microcanonical peak 
positions is approximately 90%. This value is close to 



the estimate provided by Eq. (171, which predicts a shift 
of 93% (see Table 0. 



C. Integration scheme 

We now briefly outline the numerical method used to 
evolve the system according to the equations of motion in 
Eq. Q and their integration within a standard RESPA 
multiple time scale scheme. For clarity the following 
equations are shown for a single degree of freedom. How- 
ever, extension to more dimensions follows directly. 

Integrators for molecular dynamics can be generated 
Trotter factorization of the Liouville propagator.- A suit- 
able factorization to evolve the equations of motion of a 
system coupled to a colored n oise the rmostat (Eq. ([8])) 
over a time step At is given byQSMESD 



iLAt 



s At/2 e iL p At/2 g iL x At g iL p At/2 g iL p 



,At/2 



(18) 



Each factor provides an analytic operation on the state 
of the system. The operator e* LxA * provides a coordinate 
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shift through a time step At, 

x<-x + At — (19) 
m 

while e lLpAt evolves the momenta by a time increment 
At, 



P^P + Atf(x) 



(20) 



The combination of these two operations is the stan- 
dard velocity Verlet algorithm^. The outermost opera- 
tion e lL p° At provides the effect of the colored noise ther- 
mostat on the system momentum, p and evolves the addi- 
tional thermostat d egree s of freedom, s. This operation 
can be shown to be) 18 * 20 * 



p <- Ci p + yJmk b TC 2 £ 



(21) 



where £ is a vector of independent Gaussian numbers. 
Here 



P O (22) 



is a vector containing the system momentum, p, and ex- 
tended momenta s, and 



and 



d = e -( At / 2 ) r , 



C^C2 — I ~~ C^Cj 



(23) 



(24) 



Again, it is useful to note that this reduces to a stan- 
dard integrator for the white noise Langevin thermo- 
stalpa when the extended momenta are decoupled from 



the system. In order to recover C 2 a Cholesky decompo- 
sition must be performed on the expression in Eq. ( 24 1 



This operation requires that the symmetric part of the 
matrix T is positive definite.^ 

When the system consists of more than the one de- 
gree of freedom the ability of the extended momenta to 
respond to the different frequencies experienced by each 
particle requires a local (massive) coupling!^ Hence for 
a system of N particles consisting of 6N variables (po- 
sitions and momenta) the colored noise thermostat cor- 
responding to the drift matrix in Eq. ( 14 1 adds 6N ad- 



ditional variables in the form of the auxiliary momenta, 
s. This compares favorably with those required in lo- 
cal Nose Hoover schemes which add 18N v ariables if a 
typical chain of length three is chosen EMU Evolution of 
the auxiliary momenta in Eq. (21 ) is a 3 x 3 matrix mul- 



tiplication and 3N of these operations are required for 
each thermostat evolution of the N particle system. The 
local nature of the thermostat makes the operation eas- 
ily parallelizable. For all the systems considered in this 
study we found the cost of the thermostat operations to 
be small in comparison to the force calculations which 
dominate the computational cost. 

To construct a multiple time scale scheme the forces 
are partitioned into a sum of rapidly and slowly varying 
components. In the case of three components, we may 
write the total force as, 



(25) 



where f^'(x) corresponds to the slowest forces, which 
can be integrated with the largest time step, and f^(x) 
the fastest ones which necessitate the smallest time step 
for stable integration. With this splitting of the forces 
the Liouville operator can be factorized as, 



iLAt 



^At/2 x 



Ah 

n 



e iL ps At 2 /2 e iL<, 2 »At 2 /2 



M 3 

17 (e^ 3>A * 3 



/2 e iLxAt 3e iL ( p 3 >At 3 /2\ e iL< p VAt 2 /2 e iL ps At 2 /2 



(26) 



where the exponents of Lp perform the evolution shown 



in Eq. 20 under the force f^(x). The time steps for the 



integration of the medium and fast forces are then, 



and 



At 2 = At/M 2 



At 3 = At/{M 2 M 3 ) 



(27) 



(28) 



respectively, where the whole numbers M 2 and M3 are 
chosen to be sufficiently large so as to allow stable inte- 
gration of the system under the forces (x) and /( 3 ) (x) . 



The thermostat evolution is located in the middle loop 
so that the bath is efficiently coupled to the fast system 
motions. 



III. COMPUTATIONAL DETAILS 

The colored noise RESPA scheme introduced in the 
previous section was used to perform simulations of pure 
water and alanine dipeptide in explicit water. The water 
system consisted of 1000 molecules simulated in a 31.07 A 
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FIG. 1: The shape of the memory spectra given by Eq. ( 16 1 is 
detailed schematically. In the left panel the white noise com- 
ponent (black line) is shown alongside Lorentzians centered 
at ±a)/2 (blue and red dashed curves). The right panel shows 
the memory spectra which results from the subtraction of the 
sum of the Lorentzians from the white noise component. The 
spectra is small near zero and plateaus to the value of the 
white noise component at large u. 




100 200 300 400 500 600 700 800 
CO [ ps" 1 ] 



FIG. 2: The spectra of SPC/Fw water computed from a mi- 
crocanonical simulation (solid, black line), and simulations 
that employ the GLE - 12 fs colored noise parameter set 
(red, dashed line), and a white noise (blue, dot-dashed line) 
Langevin thermostat. The memory spectra that is obtained 
from the GLE dynamics is given by the red, dotted line. The 
strength of the white noise bath corresponds to the uj — > oo 
limit of this profile (7 = 83.3 ps -1 ). 



box. The alanine dipeptide was described by the OPLS- 
AA force fiel d 32 * 33 * in a 20 A box containing 252 water 
molecules. The SPC/Fw potential was used to model 
the interactions between water molecules. All bonded 
terms (stretches, bends and torsions) were treated as flex- 
ible. 

The forces were partitioned into three levels as in Eq. 
[25] The bonded terms were updated every 0.5 fs and 
the non bonded interactions below 9 A were updated ev- 
ery 2 fs. The long range electrostatics, which dominate 
the computational cost of the simulation, were updated 
at the outer time step which was varied between 2 and 
20 fs as discussed below. The electrostatics were parti- 



tioned into short and long range parts using the RESPA2 
schemeP 3 -^ The function that switches between regions 
is a quintic spline which acts over a 4 A regionP^ The cut 
off between the short and long range electrostatics was 
chosen so as to optimize the computational performance 
within our code. Larger cutoffs have been shown to pro- 
duce better stability within the MTS formalism^ 9 so for 
maximum benefit the cutoffs used should be balanced 
with performance depending on the exact implementa- 
tion. 

Due to the cost of calculating the long range electro- 
static interactions we wish to make the outer time step 
as large as possible while still obtaining correct sampling. 
We therefore performed simulations using time steps of 
12 fs, 16 fs, and 20 fs. Resonance artifacts are extremely 
pronounced at these outer time steps so the simulation 
must be stabilized with either strong white noise damp- 
ing or optimized colored noise thermostats. The baseline 
results for dynamic and equilibrium properties were ob- 
tained from a microcanonical and a white noise Langevin 
simulation, respectively and were performed using an 
outer time step of 2 fs. For the microcanonical results 
and each combination of thermostat and time step, wa- 
ter statistics were collected over twelve independent runs 
with initial velocities sampled from the Boltzmann dis- 
tribution. For water a total of 12 ns of simulation time 
was performed, while for alanine dipeptide 400 ns was 
necessary to converge the properties reported. 



IV. RESULTS AND DISCUSSION 

For systems that exhibit large resonance artifacts, we 
now show that an appropriately designed colored noise 
thermostat is capable of yielding accurate sampling while 
simultaneously minimizing the thermostat's impact on 
diffusional and orientational dynamics. In order to test 
this scheme, we perform simulations on flexible water and 
a fully flexible simulation of aqueous alanine dipeptide. 
The broad spectrum of frequencies in aqueous systems 
range from fast intramolecular stretches to slow diffu- 
sional modes (see Figure[2]). The strong coupling between 
the modes makes this a challenging example to test our 
approach. Hence, for flexible water simulations the reso- 
nance barrier occurs at an outer time step of « 3 fs when 
the microcanonical r-RESPA algorithm is employedP 

In this work, we design colored noise thermostats that 
are capable of stabilizing resonance artifacts for outer 
time steps of 12, 16, and 20 fs while ensuring that the 
error in the energies is within 0.5% of the baseline results. 
The parameters for these thermostats are given in Table 
|TJ In order to provide comparison as to the effectiveness of 
our scheme, we also perform simulations that use a white 
noise Langevin (WNL) thermostats to yield comparable 
accuracy. These runs utilize a friction of 14.2 ps -1 , 40.0 
ps -1 , and 100.0 ps -1 in conjunction with outer times 
steps of 12 fs, 16 fs, and 20 fs, respectively. 

Additionally, recent work has shown that when rigid 
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constraints are placed on the fastest degrees of freedom, 
a weak Langevin coupling of 7 = 1 ps _1 is required to 
stabilize the simulation at an outer time step of 12 fsP 
We therefore perform simulations using a 2 f s outer time 
step with this friction. This facilitates a comparison of 
the dynamical perturbation arising from constrained dy- 
namics using white noise stabilization with that caused 
by fully flexible dynamics using our colored noise scheme. 



A. Water 

We first apply our scheme to pure flexible water. In Ta- 
ble [n] the average temperature and average bonded and 
non-bonded potential energies given for combinations of 
outer time step and method of resonance stabilization, 
either white noise Langevin (WNL) or generalized (col- 
ored noise) Langevin (GLE) thermostatting. The base- 
line result to which the static properties of all stabilized 
runs are compared utilizes an outer time step of 2 fs and 
a WNL thermostat with a friction of 7 = 1 ps _1 . It 
can be seen that the GLE runs reproduce baseline re- 
sults to within an error of 0.5% for both the bonded 
and non-bonded components of the energy. The chosen 
white noise Langevin couplings exhibit comparable over- 
all performance at each corresponding outer time step. 
Upon study of Table [TTJ it can be seen that the error 
is up to four times larger in the bonded energy as com- 
pared to the non-bonded energy for the GLE runs. In 
the case of the WNL stabilized simulations, it is up to 
ten times larger. This can be explained by the fact that 
resonance instabilities are most severe in the high fre- 
quency intramolecular modes. Therefore, as the error 
is not uniformly distributed across the system, it is of 
great utility to consider different components of the en- 
ergy when assessing resonance stabilization. Although 
the error in the temperature is slightly elevated for runs 
with an outer time step of 20 fs, it is within 0.5% for all 
other results. The errors in temperature tend to corre- 
late with those in the potential energy. However, from 
the discussion above, the potential energy can be seen to 
be a more sensitive measure of sampling accuracy. 

In Figure |3j the set of radial distributions of water of 
the GLE - 12 fs and GLE - 20 fs stabilized runs are shown 
to be in excellent agreement with the baseline calcula- 
tions. The GLE - 16 fs results for clarity are not shown 
but are fully consistent with the other results. The inset 
shows that the peak corresponding to the oxygen hydro- 
gen covalent bond is also well reproduced. This finding 
is significant since, as discussed above, the bonding en- 
ergy tends to show a larger error than the overall energy. 
This reflects the fact that these components vary rapidly, 
thereby inducing a greater sensitivity to small deviations 
in position. Overall, the results of Figure [3] underline the 
fact that the free energy surface is accurately reproduced 
by the GLE resonance stabilized simulations. 

Although both white and colored noise may be utilized 
to stabilize resonance artifacts, these two techniques sig- 



nificantly differ in the degree to which they disturb dy- 
namical properties. In Table III we report the diffusion 



constant and the relaxation time of the first order molec- 
ular dipole orientational correlation function for the reso- 
nance stabilized simulations in comparison to the micro- 
canonical baseline results. The colored noise stabilized 
dynamics with an outer time step of 12 fs differ by 4% 
from the NVE results as compared to the difference of a 
factor of 2.5 that is obtained from the WNL - 12 fs run. 
This finding is consistent with the selectivity of a colored 
noise thermostat which only couples weakly to the low 
frequency modes and therefore minimally perturbs dy- 
namical quantities that largely depend upon these slower 
modes (see Fig. [2]). It can be seen from Table III that 
the impact of the GLE - 12 fs profile on the dynamics 
is as good as and in some cases better than that of the 
weak white noise friction (7=1 ps" 1 ) which is neces- 
sary to stabilize a simulation that employs constraints at 
a 12 fs outer time stepP Therefore the present method 
of targeted damping of fast modes compares very well to 
schemes where such modes are frozen. 

The diffusion constant is extracted from the mean 
square displacement, which aside from its physical mean- 
ing, is also an indicator of the rate at which the simula- 
tion samples configuration space. Therefore, in addition 
to producing highly distorted dynamics, the strength of 
white noise coupling significantly slows down the sam- 
pling rate, largely counterbalancing any efficiency gained 
by utilizing a larger outer time step™ This drawback 
is alleviated by the use of the colored noise thermostat, 



although it can be clearly seen from Table III that the dy- 
namical perturbation engendered by our scheme increases 
with outer time step as a "stronger" colored noise profile 
is necessary to stabilize the system. This arises natu- 
rally as a greater number of modes begin to resonate at 
larger outer time steps, thereby requiring strong coupling 
to the bath across lower frequency modes to damp out 
such artifacts. In this manner, increased disturbance of 
the dynamics is necessitated. The dynamics produced 
by the GLE - 16 fs and GLE - 20 fs runs differ from 
the baseline results by approximately 20% and 35%, re- 
spectively. A comparison of the ratios of the GLE to 
the microcanonical results in Table ITTT1 confirms that this 
degree of distortion can be reasonably estimated from 
Eq. (17). However, if only equilibrium properties are de- 



sired and the computational speed-up gained by increas- 
ing the outer time step outweighs the decrease in the 
sampling rate, these choices may still offer advantages. 



B. Alanine 

Multiple time scale techniques are often employed to 
simulate biomolccular systems. In this section, we gauge 
the applicability of our scheme to a fully flexible simula- 
tion of alanine dipeptide in explicit water solvent. The 
use of this model system allows us to readily test the 
impact of the chosen colored noise profiles on conforma- 



thermostat 


outer time step (fs) 


T (K) 


Vnb/N (kcal/mol) 


Vb/N (kcal/mol) 


WNL 


2 


300.0 


-3.872 


0.5227 


GLE 


12 


301.5 


-3.869 


0.5245 


GLE 


16 


301.4 


-3.869 


0.5247 


GLE 


20 


302.0 


-3.872 


0.5236 


WNL 


12 


301.7 


-3.875 


0.5278 


WNL 


16 


301.7 


-3.874 


0.5285 


WNL 


20 


302.3 


-3.875 


0.5246 



TABLE II: Average temperatures and bonded and non-bonded potential energies per atom of the flexible water system for 
different combinations of thermostat and outer time step. The figures are reported to a precision within the error of the 
calculation. 



thermostat 


outer time step (fs) 


D (A7pb) 


Tdipolc (ps) 


NONE 


2 


0.25 


4.5 


WNL 


2 


0.22 


4.8 


GLE 


12 


0.24 


4.7 


GLE 


16 


0.20 


5.6 


GLE 


20 


0.16 


6.7 


WNL 


12 


0.095 


8.0 


WNL 


16 


0.044 


14 


WNL 


20 


0.019 


26 



TABLE III: The dynamic quantities computed in the flexi- 
ble water system. The diffusion constant and the first order 
molecular dipole relaxation time are given above. Figures are 
reported to statistical accuracy. 



tional sampling and dynamics. The colored noise profiles 
utilized are the same as those presented in Section IV A 



and here we show their applicability to more general sys- 
tems. 



In Table IV the energies of the alanine system for the 
chosen combinations of thermostats and outer time steps 
are presented. As in the case of pure water, the col- 
ored noise profiles reproduce the average bonded and 
non-bonded potential energies of the baseline (WNL - 
2 fs) result to within 0.5%. The bonded energies include 
the stretches, bends, and torsions of alanine dipeptide in 
addition to the intramolecular interactions of water. As 



also noted in Section IV A the bonded energies, which 
contain high frequency modes and most strongly exhibit 
the resonance phenomena, possess a larger overall error 
than the non-bonded energies. This behavior underlines 
why our scheme of targeting fast motions for damping 
is successful. Additionally the total and alanine dipep- 
tide temperature are reported, and can also be seen to 
be within 0.5% of the baseline results. 

The conformational space of alanine dipeptide is typi- 
cally characterized as a function of the two dihedral an- 
gles ip and <^pm!U The Ramachandran plots of the base- 
line and GLE - 12 fs result are given in Figure [4] It can be 
seen that they are in excellent agreement. The resultant 
free energy exhibits a minimum in both the a-helical and 
the extended Pn region. These two regions are labeled 
in Figure [4j 

The mean first passage time is estimated from the sur- 



- GLE - 


12fs 


- - GLE ■ 


20 fs 


... WNL 


-2fs 




2 3 
distance (A) 

FIG. 3: The atom-atom radial distribution functions of wa- 
ter computed from the baseline MTS white noise Langevin 
simulation with an an outer time step of 2 fs (black circles), 
is plotted against the colored noise thermostatted MTS pre- 
scription utilizing an outer time step of 12 fs (red solid line) 
and 20 fs (blue dashed line). The inset depicts the first peak 
of gon(r) corresponding to the OH covalent bond. 



vival probability of the transition between the an and the 
Pu regions. The survival probability, S(t), for a state in 
the an region to cross to the Pu region may be defined 
in terms of an average in conformation space over tra- 
jectories that reside in the cxr region at time t = and 
outside the Pu region up to time, t and is given by, 



S aR ^p u (t) = (l-gp u (t)) hr 



«<»)= 



(29) 



where h aR is unity inside the an region and zero outside, 
and gp u (t) is defined to be unity if the particle has passed 
into region Pu at any time between and t and is zero 
otherwise. The ip range of the a-helical region is defined 
as -40° < V < 10° and the extended Pu region as 110° < 
tp < 180°. The 4> range for both regions is set to be 
-110° <<j}< -60°. 

The survival probabilities for the an — > Pu and Pu — > 



an transitions are plotted in Figure (5 
passage times are presented in Table |V 



The mean first 
The equilibrium 



constant of the an ^ Pn process is related to the ratio 
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of the Pn — » an mean passage time to the an Pn 
value and is 1.8 — 1.9 in all runs. Upon study of the 
mean first passage times, it can be seen that the col- 
ored noise thermostated result with an outer time step 
of 12 fs is perturbed by » 10% with respect to the mi- 
crocanonical result, and is again comparable with the 
results of the system when weakly coupled to a white 
noise bath (7 = 1 ps _1 ). The colored noise thermostat 
causes a slightly larger perturbation on the mean passage 
time when compared to that exhibited in the diffusion of 
water (see Sec. IV A I. This likely arises from the depen- 
dence of this property on torsional motions that lie in 
a higher frequency range than diffusive modes, and are 
therefore more strongly coupled to the bath. The large 
white noise friction that is necessary to damp out reso- 
nance instabilities at an outer time step of 12 fs again 
increases passage times by a factor of « 2.5. The impact 
of large WNL couplings on the conformational dynam- 
ics has been noted in previous workl^SI As in the case of 
pure water, the degree of perturbation induced by the 
GLE thermostat increases as a greater degree of energy 
stabilization is necessary at larger outer time steps. The 
similarity of the results obtained in this case and pure 
water underlines the fact the thermostatting scheme pre- 
sented in this work is transferable between flexible water 
and typical solvated biomolecules. 



CONCLUSIONS 



WNL - 2 fs 




-150 -120 -90 -60 

GLE - 12 fs 




In this paper, we have shown that by coupling a stan- 
dard multiple time scale scheme to a colored noise bath 
large outer time steps can be employed while obtaining 
accurate sampling and minimally perturbing the dynam- 
ics. The scheme was illustrated with applications to flex- 
ible water and a fully flexible simulation of aqueous ala- 
nine dipeptide. Our results suggest that, when combined 
with our scheme, a 12 fs outer time step seems a good 
compromise between size of outer time step and mini- 
mizing the impact on the dynamics. With this outer 
time step the damping needed to obtain potential en- 
ergies within 0.5% of the benchmark values decreased 
diffusion by only 4%. In contrast, stabilization of the 
simulation at this outer time step using white noise de- 
creases the diffusion constant by over 2.5 times. Even 
more promising, the dynamical perturbation due to our 
colored noise scheme compared favorably with the weak 
white noise damping (7 = 1 ps _1 ) required to stabilize 
a multiple time scale simulation in which all bonds to 
hydrogen are constrainedP This suggests that, in con- 
junction with multiple time scale algorithms, this method 
may be utilized in lieu of constraints. 

In our serial code using a 12 fs outer time step yielded 
a computational speed-up of 4.3 times compared to our 
baseline MTS calculations using a 2 fs outer time step 
and a 17 times speed up compared to flexible simulations 
without a multiple time scale scheme that employ a 0.5 fs 
time step. Even further increases may be achieved when 



FIG. 4: The free energy as a function of dihedral angles (j> 
and ip. The baseline WNL - 2fs (top panel) and GLE - 12fs 
(bottom panel) results are shown. The free energy in confor- 
mation space decreases as the color varies from white to red. 
Isolines represent increments of 0.5 kcal/mol. The an and Pn 
regions are labeled in each panel. 



the computational architecture or parallelization scheme 
makes long range electrostatics comparatively costly to 
evaluate. It is therefore clear that this scheme provides 
an efficient approach to perform equilibration or config- 
urational sampling. Although our approach perturbs the 
dynamics we note that the alanine dipeptide molecule 
used in our benchmark simulation was specifically chosen 
to allow us to view many transitions and hence converge 
the dynamical properties within small bounds. In typical, 
large-scale biological system where long time-scale pro- 
cesses are of interest this would not be the case and hence 
the small change in the dynamics due to this method will 
likely be dwarfed by the statistical errors. In this case the 
ability to generate longer trajectories using this approach 
facilitates a decrease in the statistical errors bars. 

The similarities in the spectra of many aqueous and 
biological systems suggest that our noise profiles should 
provide good out-of-the box performance for other fully 
flexible systems. However, the flexibility of the colored 
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thermostat 


outer time step (fs) 


T(K) 


Tala (K) 


Vnb/N (kcal/mol) 


Vh/N (kcal/mol) 


WNL 


2 


300.0 


299.7 


-3.828 


0.5297 


GLE 


12 


301.0 


300.1 


-3.825 


0.5312 


GLE 


16 


301.1 


299.8 


-3.823 


0.5314 


GLE 


20 


301.7 


299.8 


-3.826 


0.5304 


WNL 


12 


301.3 


300.0 


-3.829 


0.5338 


WNL 


16 


301.4 


299.8 


-3.828 


0.5344 


WNL 


20 


302.0 


300.0 


-3.828 


0.5315 



TABLE IV: The total and alanine dipeptide temperatures, as well as the bonded and non-bonded energy per atom are given for 
each combination of outer time step and thermostatting scheme utilized for alanine dipeptide in explicit solvent. The figures 
are reported to statistical accuracy. 



thermostat 


outer time step (fs) 


f(a R -> Pn) (ps) 


r(P„ 


-> a R ) (ps) 


NONE 


2 


79 




142 


WNL 


2 


85 




163 


GLE 


12 


87 




156 


GLE 


16 


105 




193 


GLE 


20 


121 




234 


WNL 


12 


190 




340 


WNL 


16 


435 




790 



TABLE V: Table of mean first passage times, r, from the a R to Pn and Pn to a R regions are presented for various combinations 
of outer time steps and thermostatting schemes. The average error for runs where the dynamics are not strongly overdamped 
is ~ 3 ps. Results are not given for the WNL - 20 fs run due to the fact that too few transitions occurred in the course of the 
simulation to generate a reliable estimate of f . 




1000 



800 1200 
time (ps) 

FIG. 5: The survival probabilities for selected runs of the ala- 
nine dipeptide system are shown for the transition from the 
a r to Pn region (top panel) and for the reverse process (bot- 
tom panel). Results are given for the baseline microcanonical 
(solid black curve) and white noise Langevin (red dot-dashed 
curve) runs as well as the simulations that are resonance sta- 
bilized with an outer time step of 12 fs utilizing colored (blue 
dashed curve) and white (green dot-dashed curve) noise. 



noise approach affords such a large degree of versatil- 
ity that further tuning may yield improved results. Our 
matrices may also provide a basis for other applications 
where energy must be rapidly dissipated to ensure ac- 
curate sampling such as in the case of fast-deposition 



met adynamics.^ The methods outlined in Sec. II and 
Appendix [X] provide a starting point for such develop- 
ments. 
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Appendix A: Further applications of simple drift 
matrices 



It is often convenient to consider drift matrices whose 
elements are transparently related to the time depen- 
dence of the memory kernel (i.e. the eigenvalues of sub- 
matrix r ss ). This facilitates the use of readily tunable 
colored noise profiles that are based on a small set of pa- 
rameters. In addition to the form presented in this work 
(Eq. (14|), two other of such matrices have appeared in 
the literature. The drift matrix for a simple exponential 
noise is given 



(Al) 
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where the memory function and spectra are: 
K A (t) = 7oae -l'l a 



K a (oj) = 2 7o a 



a 



a 2 + uj a 



(A2) 
(A3) 



Whereas the following drift matrix, 



IT9I 



V 



la(a 2 +b 2 ) 
2a 

7o(a 2 +& 2 ) 
2a 



, ho(a 2 +b 2 ) l la (a 2 +b 2 ) \ 

V 2a V 2a 

a b 

-b a J 



corresponds to a term of damped oscillatory noise whose spectra is centered at ±6: 

a 2 + b 2 



K B (t) 

k B (u) 



7o- 



7o- 



-\t\a 



a 

i 2 + b 2 
a 



cos(bt) 



a 2 + {u;-bf a 2 + (uj + bf 



(A4) 



(A5) 
(A6) 



In the above equations, the parameter, a, corresponds 
to the real part and b to the magnitude of the imaginary 
part of the eigenvalue(s) of submatrix, T ss . As in Section 
PH the parameter 70 is proportional to the value of the 
memory spectra at u> = 0, such that K(0) — 270. 

Drift matrices that correspond to stable dynamics may 
be combined such that the resultant memory function is 
a weighted sum of the corresponding memory functions 
of each component.^ For example, the drift matrix that 
corresponds to the sum of the noise profiles of whose drift 
matrices are given by X T and 2 T is given by, 



(1+2) r = 



dpp + 2 Ipp) 

Isp 
Isp 




where components are expressed in the notation of 
Eq. (10 1 . In this way, forms such as those presented 
here may be utilized as "building blocks" for more flexi- 
ble colored noise profiles. 

Making use of this machinery, it is possible to design a 
form that couples strongly to certain modes and weakly 
to others. Eq. (A7) may be applied to add together noise 



Appendix B: Dynamical properties of a GLE 
dynamics 

In order to estimate the perturbation on the dynam- 
ical properties of a system caused by coupling to a col- 
ored noise thermostat one can study the behavior of a 
one-dimensional harmonic model. In this limit, micro- 
canonical dynamical results in a 5-function power spec- 
trum peaked at the oscillator's characteristic frequency 
ujq . The presence of the noise will modify the line shape 
of the peak, shifting its center and broadening it. A 
measure of these two effects can be used to estimate the 
magnitude of the disturbance. 

To achieve this, one must consider the matrices 
(A7) Y qp (u>q) and B op (uig) which describe the dynamics in the 
full x = (q,p,s) spaced One can then find the station- 
ary covariance matrix C qp by solving T qp C qp + C qp T^ p = 
BgpB^l and compute the first order correlation matrix 
/x T (t) x (0)) and its Fourier transform, 



profiles of the damped oscillatory (Eq. (A6)) or the "can- 
celing" (Eq. (16)) type. It must be noted that this for- 



malism has drawbacks when compared to the profiles uti- 
lized in the main text. Namely, the dimensionality of the 
resultant drift matrix grows with the number of targeted 
modes and that the greater sensitivity to the details of 
the spectra implies less robust performance across dif- 
ferent systems. However, this sensitivity also facilitates 
fine control of the colored noise profile for very specific 
applications. 



Cij (uj) 



F 2 +0J 2^"P 
qp ' ^ 



(CfflOii (Cqp) 



qpjjj 



-1/2 



.(Bl) 



The position and value at the maximum of the C pp (ui) 
term can then be used to characterize the deformed peak. 
When T qp is built out of the 3x3 canceling noise 



matrix in Eq. ( 14 ) and we set 70 = so as to consider 
the case of small disturbance on the low-frequency modes 
we obtain, 



■ ip 



1 


-1 





° , \ 


UJ 2 


7 


v3\/7w 









Vow 


uj 


V 


-Vt^/v/3 


— Co 


/ 



(B2) 
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Solving for C pp (w) we obtain 



3 7 w 6 

7 2 w 4 (3a; 2 + w 2 ) + (w 2 - w 2 ) (2w 2 + w 2 ) + 3 (w 2 - w 2 ) 2 (cj 4 + w 2 cj 2 + w 4 ) 



c ppH = o ~. 7z „ — -o, — /= - o ; o ^tti n — rr^ ns2 ; , — ; ( B3 ) 



<~\4 



We then take ui to be much larger than all other fre- and 
quencies entering our problem, and write them as a ratio 
with respect to Q, i.e. luq <— ujo/Q, 7 <— j/cj, to <— lo/uj. 

One then finds the relevant extremal point, ojp, and an . _ ^(uiq/Q) 

estimate of the peak width as Ap — 1/ [~kC pp (ojp)]. This ' / ■y/Cj\ 

leads to expressions which can then be expanded in pow- n \ Vs ) 

ers of wq, eventually yielding 



lop = UQ + O ({u Q /uf) (B4) 



iL 

V3 
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